### Produced Figure 3 in the paper - by pooling data from study #2 and #3.
### Also produced the alternative version in the Supplmental Information (with non attentive respondents)


### Brazil (read and combine data from the two studies)
setwd("~/Dropbox/Data/Paper-GDN/Analysis-Paper-01/Replic") #change to local folder
savedir <- getwd()

load("S02-dataBrazil.RData")
d1 <- subset(dbrazil,select=c(OUTsupport,TO,Tcond,CCTbeneficiary,attention,gender,etid))

load("S03-dataBrazil.RData")
d2 <- dbrazil
d2 <- subset(d2,select=c(OUTsupport,TO,Tcond,CCTbeneficiary,attention,gender,etid))
d <- rbind(d1,d2) 

### Examine Randomization
table(d$TO,d$Tcond)

##### Basic Barplots ###############################

robse <- function(x,sig=0.95){
	require(sandwich)
	require(car)
	robSE <- sqrt(diag(vcovHC(x,type="HC2")))
	terms <- length(coef(x))
	robP <- rep(NA,terms)
	for(jj in 1:terms){
			robP[jj] <- linearHypothesis(x,diag(terms)[jj,],vcov=vcovHC(x,type="HC2"))$Pr[2]}
	upr <- coef(x) + qnorm(sig) * robSE
	lwr <- coef(x) - qnorm(sig) * robSE
	out <- cbind(robse=robSE,upr=upr,lwr=lwr,robpval=robP)
	return(out)
}


reg02 <-lm(OUTsupport~TO*Tcond,data=subset(d,CCTbeneficiary==F&attention==T))

other.order <- c("Racial","Control","Regional")
cond.order <- c("Unconditional","Conditional")

Puncond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("_Uncond",3),
			TO=other.order),
			interval="confidence"))
Pcond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("Cond",3),
			TO=other.order),
			interval="confidence"))

to.plot <- rbind(Conditional=Pcond$fit,
				Unconditional=Puncond$fit)
colnames(to.plot)<-other.order
to.plot <- to.plot[cond.order,]


# standardized effect
effectssd <- c(racial=round((coef(reg02)[4]+coef(reg02)[5])/sd(reg02$model$OUTsupport,na.rm=T),2),
			control=round((coef(reg02)[4])/sd(reg02$model$OUTsupport,na.rm=T),2),
			regional=round((coef(reg02)[4]+coef(reg02)[6])/sd(reg02$model$OUTsupport,na.rm=T),2) )#effects in SD
	
# p-value of conditional relative unconditional		
ps <- round(c(racial=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Racial"))$p.value,
		control=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Control"))$p.value,
		regional=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Regional"))$p.value),3)
 
# p-value of premium in control relative to racial and regional
pps <-  round(robse(reg02)[c(5,4,6),"robpval"],3)
pps[2] <- NA


## Produce data for the combined effects plots (produced at the end of the code)
reg02a <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Control"))
reg02b <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Racial"))
reg02c <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Regional"))
effs <- data.frame(rbind(Control=summary(reg02a)$coef[2,1:2],
	  				Racial=summary(reg02b)$coef[2,1:2],
	  				Regional=summary(reg02c)$coef[2,1:2]))
effs$low <- effs$Estimate + qnorm(0.05)* effs$Std..Error
effs$high <- effs$Estimate + qnorm(0.95)* effs$Std..Error
effs.brazil <- effs
 
## Produce data for the combined effects plots (produced at the end of the code), this time with non attentive
## This table is included in the SI of the JDS paper
reg02aall <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&TO=="Control"))
reg02ball <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&TO=="Racial"))
reg02call <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&TO=="Regional"))
effsall <- data.frame(rbind(Control=summary(reg02aall)$coef[2,1:2],
	  				Racial=summary(reg02ball)$coef[2,1:2],
	  				Regional=summary(reg02call)$coef[2,1:2]))
effsall$low <- effsall$Estimate + qnorm(0.05)* effsall$Std..Error
effsall$high <- effsall$Estimate + qnorm(0.95)* effsall$Std..Error
effsall.brazil <- effsall



######### TURKEY

load("S02-dataTurkey.RData")
d1 <- subset(dturk,select=c(OUTsupport,TO,Tcond,CCTbeneficiary,attention,gender,etid))

load("S03-dataTurkey.RData")
d2 <- dturk
d2 <- subset(d2,select=c(OUTsupport,TO,Tcond,CCTbeneficiary,attention,gender,etid))
d <- rbind(d1,d2) 

### Examine Randomization
table(d$TO,d$Tcond)

##### Basic Barplots ###############################

reg02 <-lm(OUTsupport~TO*Tcond,data=subset(d,CCTbeneficiary==F&attention==T))

other.order <- c("Racial","Control","Regional")
cond.order <- c("Unconditional","Conditional")

Puncond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("_Uncond",3),
			TO=other.order),
			interval="confidence"))
Pcond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("Cond",3),
			TO=other.order),
			interval="confidence"))

to.plot <- rbind(Conditional=Pcond$fit,
				Unconditional=Puncond$fit)
colnames(to.plot)<-other.order
to.plot <- to.plot[cond.order,]


# standardized effect
effectssd <- c(racial=round((coef(reg02)[4]+coef(reg02)[5])/sd(reg02$model$OUTsupport,na.rm=T),2),
			control=round((coef(reg02)[4])/sd(reg02$model$OUTsupport,na.rm=T),2),
			regional=round((coef(reg02)[4]+coef(reg02)[6])/sd(reg02$model$OUTsupport,na.rm=T),2) )#effects in SD
	
# p-value of conditional relative unconditional		
ps <- round(c(racial=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Racial"))$p.value,
		control=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Control"))$p.value,
		regional=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Regional"))$p.value),3)
 
# p-value of premium in control relative to racial and regional
pps <-  round(robse(reg02)[c(5,4,6),"robpval"],3)
pps[2] <- NA


## Produce data for the combined efc(2,3,1)fects plots (produced at the end of the code)
reg02a <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Control"))
reg02b <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Racial"))
reg02c <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Regional"))
effs <- data.frame(rbind(Control=summary(reg02a)$coef[2,1:2],
	  				Racial=summary(reg02b)$coef[2,1:2],
	  				Regional=summary(reg02c)$coef[2,1:2]))
effs$low <- effs$Estimate + qnorm(0.05)* effs$Std..Error
effs$high <- effs$Estimate + qnorm(0.95)* effs$Std..Error
effs.turkey <- effs

## Produce data for the combined effects plots (produced at the end of the code), this time with non attentive
## This table is included in the SI of the JSD paper
reg02aall <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&TO=="Control"))
reg02ball <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&TO=="Racial"))
reg02call <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&TO=="Regional"))
effsall <- data.frame(rbind(Control=summary(reg02aall)$coef[2,1:2],
	  				Racial=summary(reg02ball)$coef[2,1:2],
	  				Regional=summary(reg02call)$coef[2,1:2]))
effsall$low <- effsall$Estimate + qnorm(0.05)* effsall$Std..Error
effsall$high <- effsall$Estimate + qnorm(0.95)* effsall$Std..Error
effsall.turkey <- effsall


### The pooled plot - FIgure 3 in the paper #####
pdf(file=paste(savedir,"/fig-P01S0102-pooledeffects.pdf",sep=""))
plot(c(-.2,1),c(0.5,3.5),type="n",bty="n",ylab="",xlab="",yaxt="n")
#polygon(x=c(-.2,1,1,-.2),y=c(1.5,1.5,2.5,2.5),border=NA,col=gray(0.7))
abline(v=0,lty=3)
segments(x0=effs.brazil$low,x1=effs.brazil$high,y0=c(2,3,1)+0.15,y1=c(2,3,1)+0.15)
segments(x0=effs.turkey$low,x1=effs.turkey$high,y0=c(2,3,1)-0.15,y1=c(2,3,1)-0.15)
points(effs.brazil$Estimate,c(2,3,1)+0.15,pch=21,bg="white")
points(effs.turkey$Estimate,c(2,3,1)-0.15,pch=24,bg="black")
legend(x="topright",bty="n",legend=c("Brazil","Turkey"),pch=c(21,24),pt.bg=c("white","black"),inset=0.015)
axis(side=2,at=c(1,2,3),labels=c("Regional","Control","Racial"),las=2,tick=F,hadj=0.5)
mtext(side=1,line=3,"Standardized Effects")
dev.off()


### The pooled plot including non-attentive respondents in the Supplemental Information of the JDS paper #####
pdf(file=paste(savedir,"/fig-P01S0102-pooledeffects-all.pdf",sep=""))
plot(c(-.2,1),c(0.5,3.5),type="n",bty="n",ylab="",xlab="",yaxt="n")
#polygon(x=c(-.2,1,1,-.2),y=c(1.5,1.5,2.5,2.5),border=NA,col=gray(0.7))
abline(v=0,lty=3)
segments(x0=effsall.brazil$low,x1=effsall.brazil$high,y0=c(2,3,1)+0.15,y1=c(2,3,1)+0.15)
segments(x0=effsall.turkey$low,x1=effsall.turkey$high,y0=c(2,3,1)-0.15,y1=c(2,3,1)-0.15)
points(effsall.brazil$Estimate,c(2,3,1)+0.15,pch=21,bg="white")
points(effsall.turkey$Estimate,c(2,3,1)-0.15,pch=24,bg="black")
legend(x="topright",bty="n",legend=c("Brazil","Turkey"),pch=c(21,24),pt.bg=c("white","black"),inset=0.015)
axis(side=2,at=c(1,2,3),labels=c("Regional","Control","Racial"),las=2,tick=F,hadj=0.5)
mtext(side=1,line=3,"Standardized Effects")
dev.off()




#### Equivalent data for Chile and Uruguay, not included in paper 1 (requested by GDN) #####



#### CHILE 
setwd("~/Dropbox/Data/Paper-GDN/Analysis-Paper-01")
savedir <- "~/Dropbox/LatexFiles/Paper-GDN/FiguresPaper01"
load("~/Dropbox/data/Paper-GDN/dataChile.RData")
d1 <- subset(dchile,Tcond!="Cond")

d1$Tcond <- ifelse(d1$Tcond=="CondChild","Cond",d1$Tcond)
d1 <- subset(d1,select=c(OUTsupport,TO,Tcond,CCTbeneficiary,attention))

load("~/Dropbox/data/Paper-GDN/dataChileDirect.RData")
d2 <- dchile
d2 <- subset(d2,select=c(OUTsupport,TO,Tcond,CCTbeneficiary,attention))
d <- rbind(d1,d2) 

### Examine Randomization
table(d$TO,d$Tcond)

##### Basic Barplots ###############################
sel <- 3

### Interaction Effects ####
robse <- function(x,sig=0.95){
	require(sandwich)
	require(car)
	robSE <- sqrt(diag(vcovHC(x,type="HC2")))
	terms <- length(coef(x))
	robP <- rep(NA,terms)
	for(jj in 1:terms){
			robP[jj] <- linearHypothesis(x,diag(terms)[jj,],vcov=vcovHC(x,type="HC2"))$Pr[2]}
	upr <- coef(x) + qnorm(sig) * robSE
	lwr <- coef(x) - qnorm(sig) * robSE
	out <- cbind(robse=robSE,upr=upr,lwr=lwr,robpval=robP)
	return(out)
}


if(sel==1){reg02 <- lm(OUTsupport~TO*Tcond,data=subset(d,CCTbeneficiary==F))}
if(sel==2){reg02 <- lm(OUTsupport~TO*Tcond,data=subset(d,attention==T))}
if(sel==3){reg02 <-lm(OUTsupport~TO*Tcond,data=subset(d,CCTbeneficiary==F&attention==T))}

other.order <- c("Racial","Control","Regional")
cond.order <- c("Unconditional","Conditional")

Puncond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("_Uncond",3),
			TO=other.order),
			interval="confidence"))
Pcond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("Cond",3),
			TO=other.order),
			interval="confidence"))

to.plot <- rbind(Conditional=Pcond$fit,
				Unconditional=Puncond$fit)
colnames(to.plot)<-other.order
to.plot <- to.plot[cond.order,]


# standardized effect
effectssd <- c(racial=round((coef(reg02)[4]+coef(reg02)[5])/sd(reg02$model$OUTsupport,na.rm=T),2),
			control=round((coef(reg02)[4])/sd(reg02$model$OUTsupport,na.rm=T),2),
			regional=round((coef(reg02)[4]+coef(reg02)[6])/sd(reg02$model$OUTsupport,na.rm=T),2) )#effects in SD
	
# p-value of conditional relative unconditional		
ps <- round(c(racial=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Racial"))$p.value,
		control=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Control"))$p.value,
		regional=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Regional"))$p.value),3)
 
# p-value of premium in control relative to racial and regional
pps <-  round(robse(reg02)[c(5,4,6),"robpval"],3)
pps[2] <- NA

## Produce data for the combined efc(2,3,1)fects plots (produced at the end of the code)
reg02a <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Control"))
reg02b <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Racial"))
reg02c <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Regional"))
effs <- data.frame(rbind(Control=summary(reg02a)$coef[2,1:2],
	  				Racial=summary(reg02b)$coef[2,1:2],
	  				Regional=summary(reg02c)$coef[2,1:2]))
effs$low <- effs$Estimate + qnorm(0.05)* effs$Std..Error
effs$high <- effs$Estimate + qnorm(0.95)* effs$Std..Error
effs.chile<- effs



##### Uruguay

#### CHILE 
setwd("~/Dropbox/Data/Paper-GDN/Analysis-Paper-01")
savedir <- "~/Dropbox/LatexFiles/Paper-GDN/FiguresPaper01"
load("~/Dropbox/data/Paper-GDN/dataUruguay.RData")
d1 <- subset(dchile,Tcond!="Cond")

d1$Tcond <- ifelse(d1$Tcond=="CondChild","Cond",d1$Tcond)
d1 <- subset(d1,select=c(OUTsupport,TO,Tcond,CCTbeneficiary,attention))

load("~/Dropbox/data/Paper-GDN/dataUruguayDirect.RData")
d2 <- dchile
d2 <- subset(d2,select=c(OUTsupport,TO,Tcond,CCTbeneficiary,attention))
d <- rbind(d1,d2) 

### Examine Randomization
table(d$TO,d$Tcond)

##### Basic Barplots ###############################
sel <- 3

### Interaction Effects ####
robse <- function(x,sig=0.95){
	require(sandwich)
	require(car)
	robSE <- sqrt(diag(vcovHC(x,type="HC2")))
	terms <- length(coef(x))
	robP <- rep(NA,terms)
	for(jj in 1:terms){
			robP[jj] <- linearHypothesis(x,diag(terms)[jj,],vcov=vcovHC(x,type="HC2"))$Pr[2]}
	upr <- coef(x) + qnorm(sig) * robSE
	lwr <- coef(x) - qnorm(sig) * robSE
	out <- cbind(robse=robSE,upr=upr,lwr=lwr,robpval=robP)
	return(out)
}


if(sel==1){reg02 <- lm(OUTsupport~TO*Tcond,data=subset(d,CCTbeneficiary==F))}
if(sel==2){reg02 <- lm(OUTsupport~TO*Tcond,data=subset(d,attention==T))}
if(sel==3){reg02 <-lm(OUTsupport~TO*Tcond,data=subset(d,CCTbeneficiary==F&attention==T))}

other.order <- c("Racial","Control","Regional")
cond.order <- c("Unconditional","Conditional")

Puncond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("_Uncond",3),
			TO=other.order),
			interval="confidence"))
Pcond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("Cond",3),
			TO=other.order),
			interval="confidence"))

to.plot <- rbind(Conditional=Pcond$fit,
				Unconditional=Puncond$fit)
colnames(to.plot)<-other.order
to.plot <- to.plot[cond.order,]


# standardized effect
effectssd <- c(racial=round((coef(reg02)[4]+coef(reg02)[5])/sd(reg02$model$OUTsupport,na.rm=T),2),
			control=round((coef(reg02)[4])/sd(reg02$model$OUTsupport,na.rm=T),2),
			regional=round((coef(reg02)[4]+coef(reg02)[6])/sd(reg02$model$OUTsupport,na.rm=T),2) )#effects in SD
	
# p-value of conditional relative unconditional		
ps <- round(c(racial=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Racial"))$p.value,
		control=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Control"))$p.value,
		regional=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Regional"))$p.value),3)
 
# p-value of premium in control relative to racial and regional
pps <-  round(robse(reg02)[c(5,4,6),"robpval"],3)
pps[2] <- NA

## Produce data for the combined efc(2,3,1)fects plots (produced at the end of the code)
reg02a <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Control"))
reg02b <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Racial"))
reg02c <-lm(scale(OUTsupport)~Tcond,data=subset(d,CCTbeneficiary==F&attention==T&TO=="Regional"))
effs <- data.frame(rbind(Control=summary(reg02a)$coef[2,1:2],
	  				Racial=summary(reg02b)$coef[2,1:2],
	  				Regional=summary(reg02c)$coef[2,1:2]))
effs$low <- effs$Estimate + qnorm(0.05)* effs$Std..Error
effs$high <- effs$Estimate + qnorm(0.95)* effs$Std..Error
effs.uruguay<- effs

pdf(file=paste(savedir,"/fig-P01S0102h-pooledeffects.pdf",sep=""))
plot(c(-.2,1),c(0.5,3.5),type="n",bty="n",ylab="",xlab="",yaxt="n")
#polygon(x=c(-.2,1,1,-.2),y=c(1.5,1.5,2.5,2.5),border=NA,col=gray(0.7))
abline(v=0,lty=3)
segments(x0=effs.chile$low,x1=effs.chile$high,y0=c(2,3,1)+0.15,y1=c(2,3,1)+0.15)
segments(x0=effs.uruguay$low,x1=effs.uruguay$high,y0=c(2,3,1)-0.15,y1=c(2,3,1)-0.15)
points(effs.chile$Estimate,c(2,3,1)+0.15,pch=21,bg="white")
points(effs.uruguay$Estimate,c(2,3,1)-0.15,pch=24,bg="black")
legend(x="topright",bty="n",legend=c("Chile","Uruguay"),pch=c(21,24),pt.bg=c("white","black"),inset=0.015)
axis(side=2,at=c(1,2,3),labels=c("Regional","Control","Racial"),las=2,tick=F,hadj=0.5)
mtext(side=1,line=3,"Standardized Effects")
dev.off()




#Save data for GDN
write.csv(effs.brazil,file="~/Dropbox/LatexFiles/Paper-GDN/GDN-MSWord Version/data-othernessbrazil.csv")
write.csv(effs.turkey,file="~/Dropbox/LatexFiles/Paper-GDN/GDN-MSWord Version/data-othernessturkey.csv")
write.csv(effs.chile,file="~/Dropbox/LatexFiles/Paper-GDN/GDN-MSWord Version/data-othernesschile.csv")
write.csv(effs.uruguay,file="~/Dropbox/LatexFiles/Paper-GDN/GDN-MSWord Version/data-othernessuruguay.csv")



